Astronomy & Astrophysics manuscript no. 3DWavelets 


©ESO 2012 


January 23, 2012 





Spherical 3D Isotropic Wavelets 

Francois Lanusse 1 , Anais Rassat 2 ' 1 , and Jean-Luc Starck 1 

1 Laboratoire AIM, UMR CEA-CNRS-Paris 7, Irfu, SEDI-SAP, Service d'Astrophysique, CEA Saclay, F-91191 GIF-Sur-YVETTE 
CEDEX, France. 

2 Laboratoire d'Astrophysique, Ecole Polytechnique Federate de Lausanne (EPFL), Observatoire de Sauverny, CH-1290, Versoix, 
Switzerland. 

January 23, 2012 

ABSTRACT 

Context. Future cosmological surveys will provide 3D large scale structure maps with large sky coverage, for which a 3D Spherical 
Fourier-Be ssel (SFB) analysis in spherical coordinates is natural. Wavelets are particularly well-suited to the analysis and denoising 
of cosmological data, but a spherical 3D isotropic wavelet transform does not currently exist to analyse spherical 3D data. 
Aims. The aim of this paper is to present a new formalism for a spherical 3D isotropic wavelet, i.e. one based on the SFB decomposition 
of a 3D field and accompany the formalism with a public code to perform wavelet transforms. 

Methods. We describe a new 3D isotropic spherical wavelet decomposition based on the undecimated wavelet transform (UWT) 
described in Starck et al. 2006. We also present a new fast Discrete Spherical Fourier-Bessel Transform (DSFBT) based on both a 
discrete Bessel Transform and the HEALPIX angular pixelisation scheme. We test the 3D wavelet transform and as a toy-application, 
apply a denoising algorithm in wavelet space to the Virgo large box cosmological simulations and find we can successfully remove 
noise without much loss to the large scale structure. 

Results. We have described a new spherical 3D isotropic wavelet transform, ideally suited to analyse and denoise future 3D spherical 
cosmological surveys, which uses a novel Discrete Spherical Fourier-Bessel Transform. We illustrate its potential use for denoising 
using a toy model. All the algorithms presented in this paper are available for download as a public code called MRS 3D at http : 
// j starck . free . £r/mrs3d . html 

Key words. Wavelets - Spherical Fourier-Bessel, Cosmology: Large-Scale Structure of the Universe, Methods: Data Analysis, 
Methods: Statistical 



1. Introduction 

Challenges in Modem Cosmology 

The wealth of cosmological data in the last few decades (Larson 
et al. 2010; Schrabback et al. 2010; Percival et al. 2007 b) has led 
to the establishment of a standard model of cosmology, which 
describes the Universe as composed today of approximately 4% 
baryons, 22% dark matter and 74% dark energy. The main chal- 
lenges in modern cosmology are to understand the nature of both 
dark energy and dark matter, as well as the initial conditions of 
the Universe (Albrecht et al. 2006; Peacock et al. 2006). A thor- 
ough understanding of these three topics may lead to a revision 
of Einstein's theory of General Relativity and our view of the 
early Universe. 

New surveys are planned who aim to answer these important 
questions e.g. Planck for the CMB (Planck Collaboration et al. 
2011), DES (Dark Energy Survey, Annis et al. 2005), BOSS 
(Baryon Oscillation Spectroscopic Survey, Schlegel et al. 2007), 
LSST (Large Synoptic Survey Telescope, Tyson & LSST 2004) 
and Euclid (Laureijs et al. 201 1) for weak lensing and the study 
of large scale structure with galaxy surveys. 

The challenge with these upcoming large data- sets is to ex- 
tract the cosmological information in the most suitable man- 
ner in order to test the cosmological paradigm. Depending on 
the signal one wishes to extract, and/or survey geometry, differ- 
ent bases may be more or less suitable (e.g., Fourier, Spherical 
Harmonic, Configuration or Wavelet Space). Moreover, future 
surveys may be in 2D (e.g. Planck) or in 3D (e.g. galaxy or weak 
lensing surveys). Where 3D data is available, a tomographic 



analysis is possible (also known as 2D 1/2), or a full 3D analysis 
can be done. For data in spherical coordinates, this corresponds 
to a Spherical Fourier-Bessel (SFB) decomposition (Heavens & 
Taylor 1995; Fisher et al. 1995; Castro et al. 2005; Erdogdu et al. 
2006a,b; Leistedt et al. 2011; Rassat & Refregier 2011). 

Wavelet Transform on the Sphere 

Wavelets are particularly well suited to the analysis of cosmo- 
logical data (Martinez et al. 1993; Starck et al. 2006), since 
cosmological data can often be sparsely represented in wavelet 
space. 2D Wavelets have been used in many as trophy sical stud- 
ies (Starck & Murtagh 2006) for a broad range of applications 
such as denoising, deconvolution, detection, etc. CMB studies 
have motivated the development of 2D spherical wavelet decom- 
positions. Continuous wavelet transforms on the sphere (Antoine 
1999; Tenorio et al. 1999; Cayon et al. 2001; Holschneider 1996) 
have been proposed, mainly for non Gaussianity studies. In 
Starck et al. (2006), an invertible isotropic undecimated wavelet 
transform (UWT) on the sphere based on spherical harmonics 
was described, that can be also used for other applications such 
as deconvolution, component separation (Moudden et al. 2005; 
Bobin et al. 2008; Delabrouille et al. 2008), inpainting (Abrial 
et al. 2007; Abrial et al. 2008), or Poisson denoising (Schmitt 
et al. 2010). A similar wavelet construction has been published 
in (Marinucci et al. 2008; Fay & Guilloux 2008; Fay et al. 2008) 
using so-called "needlet filters", and in Wiaux et al. (2008), an 
algorithm was proposed which allows us to reconstruct an im- 
age from its steerable wavelet transform. Other multiscale trans- 
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forms on the sphere such as ridgelets and curvelets have been 
developed (Starck et al. 2006), and are well adapted to detect 
anisotropic features. Other multiscale transforms on the sphere, 
such as ridgelets and curvelets, have been developed (Starck 
et al. 2006) and are well adapted to detect anisotropic features. 
An extension of this UWT has also been developed for polarised 
CMB data in Starck et al. (2009). 

In this paper, we describe a new 3D isotropic spherical 
wavelet decomposition, which is reversible, and could therefore 
be useful for many different applications. It is based on the UWT 
proposed by Starck et al. (2006) and extended into 3D. The 3D- 
UWT proposed here can be used to analyse 3D data in spherical 
coordinates, such as a 3D galaxy or weak lensing survey with 
large (but not necessarily full) sky coverage. 

2. Spherical 3D Filtering using the Spherical 
Fourier-Bessel Transform 

2.1. Spherical Fourier-Bessel Transform 

The SFB transform of a square integrable scalar field f(r, 6, (p) 
can be defined as: 



flm(k) 



= ^ JJJ f(r,6,^)ji(kr)Y^(6^)r 2 sin(6)drd6d(p, 



(1) 

where Y™ are spherical harmonics and ji are spherical Bessel 

functions and Y represents the complex conjugate of Y. Note 
Equation 1 uses the same convention as Heavens & Taylor 
(1995) which differs slightly from that of Castro et al. (2005), 
Leistedt et al. (2011) and Rassat & Refregier (2011), as ex- 
plained below. This expression allows the expansion of a 3D 
field provided in spherical coordinates onto a set of orthogonal 
functions: 



Mkr)Y™(6,cf>) 



(2) 



From the orthogonality of the ^/^'s, the original field 
/(r, 6, cp) can be reconstructed from its SFB transform fi m (k) us- 
ing the following inversion formula: 

a/^ZZ ffim(k)k 2 ji(kr)dkY lm (6,(P). (3) 

V 71 1=0 m=-l J 



f(r,6, 



From this definition, the SFB transform can also be regarded 
as the commutative composition of two different transforms: a 
Spherical Harmonics Transform (SHT) for the angular dimen- 
sion and a Spherical Bessel Transform (SBT) for the radial di- 
mension. In this work, the following convention is adopted to 
define the SHT of a function f(0, (p) defined on the sphere: 

flm 



r»Z7T r*7T 

Jo Jo 



Yfie, <f>)f(0, 0) sin(6)d6d(p, (4a) 



1=0 m=-l 



(4b) 



Our choice for Y; mJ t allows us to give a symmetrical expres- 
sion for the SBT and its inverse: 



fi(k) = 
f(r) = 





f(r)ji(kr)r 2 dr, 



(5a) 



(5b) 



Although this symmetrical formulation for the SBT may differ 
from other works (e.g., Castro et al. 2005; Leistedt et al. 2011; 
Rassat & Refregier 2011), it will prove very convenient, espe- 
cially to obtain a discretised transform (see Section 4). 

2.2. Frequency filtering using the SFB Transform 

Spherical 3D filtering can be defined as the 3D convolution prod- 
uct of a 3D field in spherical coordinates with a 3D filter also 
provided in spherical coordinates. Using the relations presented 
in A.l and A.2, it is possible to express such a product in terms 
of SFB coefficients and to relate those coefficients to regular 
Fourier coefficients. 

To illustrate Spherical 3D filtering on a simple case, let us 
consider a simple isotropic low-pass filter u whose 3D Fourier 
transform is U(k, 0^). The 3D Fourier transform of such a fil- 
ter has a spherical symmetry and takes a simple form in spherical 
coordinates. Indeed, U(k, 6k, (pk) is a function only of k because 
of its symmetry, therefore U(k,6k,(pk) = U(k). Furthermore, 
since U (k, 6k, (pk) is constant on every sphere centred on the ori- 
gin in Fourier space, its SHT for a given k verifies V(/,m) ^ 
(0, 0), Ui m (k) = 0. As a result, the SFB transform of u has the 
following simple expression: 



iiimik) 



V -f if/ = m = 0, 







otherwise. 



(6) 



Let us now consider a 3D field / defined in spherical co- 
ordinates. The filtered field v resulting from applying u to / is 
obtained using Eq. (A. 8) simplified by the properties of g\ 



vim(k) = (f*u) lm (k) 



(i/V^^ 2 (-i) V fi>m>(k) 



Z'=0 m'=-V 



i+r 

Z c l \l m, /', m')(-i) 1 "ui» m - m > (k)6i»o6 m »o 9 



i"=\i-v\ 



= H 2 )' V(2^i(^ (/, m, /, m)(-0°fioo(*). 



(7) 



Knowing that c°(l, m, /, m) = 1/ V47T, the following expression is 
finally obtained: 



t (k) = ^27iuoo{k)fi m {k). 



(8) 



In the special case of a 3D isotropic filter, frequency filtering 
is therefore easily obtained using the SFB transform and the fil- 
tered coefficients are simply the original coefficients multiplied 
by a function of k. 

Although this filter seems overly simplistic, it will be shown 
in the following sections that such a low-pass filter is at the heart 
of the Isotropic Undecimated Spherical 3D Wavelet transform 
which makes direct use of Eq. (8). 

3. Isotropic Undecimated Spherical 3D Wavelet 
Transform 

An Isotropic Undecimated Spherical Wavelet transform defined 
on the sphere and based on the SHT was proposed in Starck et al. 
(2006). Now the aim is to transpose the ideas behind this trans- 
form to the case of data in 3D spherical coordinates. Indeed, 
the isotropic wavelet transform can be defined using only an 
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isotropic low-pass filter. In the last section the necessary rela- 
tions to apply such a filter have been obtained and the Isotropic 
Undecimated Spherical 3D Wavelet transform can now be de- 
fined. 



3.1. Wavelet decomposition 

Using the formalism introduced in the previous section, a 
Wavelet transform can be defined with the SFB transform. This 
isotropic transform is based on a scaling function fi kc (r,6 r ,(p r ) 
with cut-off frequency k c and spherical symmetry. The symmetry 
of this function is preserved in the Fourier space and therefore, 
its SFB transform verifies = whenever (/, m) ± (0, 0). 

Furthermore, due to its cut-off frequency, the scaling function 
verifies (f>Q (k) = for all k > & c . In other terms, the scaling 
function verifies: 



&°(r,e n ipr) = &°{r) = A /- 




® k w(k)k 2 Mkr)dkY°(e r ,<t> r ). 

(9) 

Using relation (8) the convolution of the original data f(r, 0, (p) 
with <S> kc becomes very simple: 



c° 



lm (k) = [® kc *f] lm (k) = ^2n¥^(k)f lm {k). (10) 

Using this scaling function, it is possible to define a sequence 
of smoother approximations c j (r, r , cp r ) of a function /(r, r , (f> r ) 
on a dyadic resolution scale. Let O 2 ' kc be a rescaled version of 
<$> kc with cut-off frequency 2~ j k c . Then c j (r, r , <f) r ) is obtained by 
convolving f(r, r , 4> r ) with 5> 2 ' kc : 



c l = 0> 2 ^* 



f, 



Applying the SFB transform to the last relation yields: 



^/2n^(k)f lm (k). 



(11) 



(12) 



This leads to the following recurrence formula : 



By applying the SFB transform to the definition of the 
wavelet coefficients and using the recurrence formula verified 
by the c j s yields: 



Vk<^, 



^ lm ^ 



"lm 



(*). (17) 



Equations (17) et (13) which define the wavelet decomposi- 
tion are in fact equivalent to convolving the resolution at a given 
scale j with a low-pass and a high-pass filter in order to obtain 
respectively the resolution and the wavelet coefficients at scale 
7+1. 

The low-pass filter, denoted here by h j can be defined for 
each scale j by : 



(k) 



ifk< 



k c 







2-i +l 

otherwise. 



and / : 



0, 



(18) 



Then the approximation at scale j + 1 is given by the convo- 
lution of scale j with W : 



1 



; Q J * 



y/27r 



-h J . 



(19) 



In the same way, a high pass filter g j can be defined on each 
scale j by: 



Slik) 



1 

10 



if k < and / = m = 0, 

if k > and / = m = 0, 
otherwise. 



(20) 



Given the definition of X F, g j can also be expressed in the simple 
form : 

sL(k) = l-hi m (k). (21) 

The wavelet coefficients at scale j are obtained by convolving 
the resolution at scale j - 1 with g j ~ l \ 



w 7 = c ] 1 * 



V2tt 



yj- 1 



(22) 



c j+1 ( 

lm 



^00 



"ImS 



(13) 



Just like with the a trous algorithm by Starck et al. (2010), 
the wavelet coefficients {w 7 } can now be defined as the difference 
between two consecutive resolutions: 



In summary, the two relations necessary to recursively define 
the wavelet transform are: 



S J m (k)ci m (k). 



(23) 
(24) 



w j+l (r, e r , <f> r ) = c J (r, 0, <j>) - c J+l (r, 0, <f>). 



(14) 



This choice for the wavelet coefficients is equivalent to the 
following definition for the wavelet function *¥ kc : 



lm 

so that : 

= m k °*f, 

1 _ vp2->i^ /; 



(k) 



lm 



w 



w 



= * f m 



(15) 



(16) 



3.2. Reconstruction 

Since the wavelet coefficients are defined as the difference be- 
tween two resolutions, the reconstruction from the wavelet de- 



composition <W - {w 



is straightforward and corre- 



sponds to the reconstruction formula of the a trous algorithm: 

j 



C°(r, 6 r , (f> r ) = C J (r, 6 r , (f> r ) + ^ °r, <frr\ 

7=1 



(25) 



However, given the redundancy of the transform, the reconstruc- 
tion is not unique. It is possible to take advantage of this re- 
dundancy to reconstruct c j from c 7+1 and w 7+1 by using a least 
squares estimate. 
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From the previous paragraph, the wavelet decomposition can 
be then recursively defined by: 



(26) 
(27) 



If these relations are respectively multiplied by h J [m (k) and g J lm (k) 
and then added together, the following expression is obtained for 
the least squares estimate of d from c J+1 and w j+1 : 



where V and g-> are defined as follows: 



h j i m (k) 



\hL(k)\ 2 + \gijk)\ 



hn 



(k) 



(28) 

(29) 
(30) 



Among the advantages of using this reconstruction formula 
instead of the raw sum over the wavelet coefficients is that there 
is no need to perform an inverse and then direct SFB transform 
to reconstruct the coefficients of the original data. Indeed, both 
the wavelet decomposition and reconstruction procedures only 
require access to SFB coefficients and there is no need to revert 
back to the direct space (although this will be required if one 
wishes to observe the power of each coefficient in real space as 
we do in Section 5 for denoising purposes). 

3.3. Choice of a scaling function 

Any function with spherical symmetry and a cut-off frequency 
k c would do as a scaling function but in this work we choose to 
use a B- spline function of order 3 to define our scaling function: 



3 / 2& 



(31) 



where 

^3 W = ^ (Ix - 2,3 - 4|, - IP + 6|,|3 - 4|x + IP + „ + 2|3) . 

(32) 

The scaling function and its corresponding wavelet function 
are plotted in SFB space for different values of j in Figure 1 . 

Other functions such as Meyer wavelets or the needlet func- 
tion (Marinucci et al. 2008) can be used as well. Needlet wavelet 
functions have a much better frequency localisation than the 
wavelet function derived from the B 3 -spline but present more 
oscillations in the direct space. To illustrate this, we show in 
Fig. 2 two different wavelet functions. Fig. 2 left shows ID 
profile of the spline (continuous line) and the needlet wavelet 
function (dotted line) at a given scale. Fig. 2 right shows the 
same function, but we have plotted the absolute value in order 
to better visualise their respective ringing. As it can be seen, for 
wavelet functions with the same main lob, the needlet wavelet 
oscillates much more than the spline wavelet. Hence, the best 
wavelet choice certainly depends on the desired applications. For 
statistical analysis, detection or restoration applications, we may 
prefer to use a wavelet which does not oscillate too much and 
with a smaller support, in this case the spline wavelet is clearly 
the correct choice. For spectral or bispectral analysis, where the 
frequency localisation is fundamental, then needlets should be 
preferred to the spline wavelet. 



4. The Discrete Spherical Fourier-Bessel Transform 
(DSFBT) 

The spherical 3D wavelet formalism derived in the previous sec- 
tion is based on a continuous SFB transform. However, a contin- 
uous transform is not practical numerically and a discrete equiv- 
alent is required. In the case of galaxy surveys, where the un- 
derlying density field is sampled at discrete points, this problem 
can be addressed by combining a boundary condition with the 
discrete nature of the survey (see for e.g. Equations 5 and 7 in 
Leistedt et al. 201 1, even though they use a slightly different SFB 
expansion as in this paper) which yields in this specific case a 
practical way to estimate discrete SFB coefficients. However this 
method cannot be used for the purpose of performing wavelet 
treatments since this requires a continuous field. 

Indeed, to perform operations on the wavelet scales (e.g. for 
denoising purposes), one first needs to recover the wavelet co- 
efficients in physical space, apply a treatment to those coeffi- 
cients (e.g. thresholding) and then convert them back into SFB 
space (as required by the wavelet transform algorithm which op- 
erates in SFB space). This last operation requires a SFB trans- 
form which takes as input fields and not discrete galaxy surveys. 

In this section, we show how a discretisation of the k spec- 
trum in combination with a Healpix type angular pixelisation, 
can be used to build a fast Discrete SFB Transform (DSFBT) 
which allows for a novel and practical two-way conversion be- 
tween the discrete SFB coefficients and discrete values of the 
field in physical space sampled on a spherical 3D grid. 

4.1. The Discrete Spherical Bessel Transform 

The discretisation of the SFB transform along the radial compo- 
nent uses the well known orthogonality property of the Spherical 
Bessel functions on the interval [0, R]. If / is a continuous func- 
tion defined on [0,R] which verifies the boundary condition 
f(R) = then the SBT defined by Eq. (5) can be expressed using 
SFB series: 



fiihn) 




f(r)ji(ki n r)r 2 dr, 



x Jo 

00 

f(r) = ^jfi(ki n )pi n ji(ki n r). 



(33) 
(34) 



77=1 



In this expression, k\ n - ^ where q\ n is the nth zero of the Bessel 
function of order / and the weights p\ n are defined as: 



V2^r 3 

jf + Mln) 



(35) 



Although this formulation provides a discretisation of the in- 
verse SBT and of the k spectrum, the direct transform is still 
continuous and another discretisation step is required. 

If a boundary condition of the same kind is applied to fi(k) so 
that fi(Ki) = 0, then by using the same result, the SFB expansion 
of fi(k) is obtained by: 



fiiTin) = 

fm = 




fi(k)ji(n n k)k dk, 



(36a) 



2 fi(rin)Kinji(ri n k), (36b) 



71=1 
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:i.«:= 



i ji 



1.7 



2.:j 



(a) Scaling function On ^ c (£) for j = 0, 1, 2 




(b) Wavelet function ¥™ f o r 7 = 0, 1, 2 



Fig. 1. Scaling function and Wavelet function for k c = 1. The y-axis represents the amplitude and the x-axis the frequency associated 
with the scaling function. The units of the amplitude must be those of the SFB transform. 



Spline wavelet (continuous) versus Needlet wavelet (dotted) ABS(Spline wavelet) (continuous) versus ABS(Needlet wavelet) (dotl 




Fig. 2. Comparaison between spline and needlet wavelet functions on the sphere. The j-axis corresponds to amplitude and the x-axis 
to position. 



where r/ n = ®r and where the weights K\ n are defined as: 



Kin 



j 2 i+Min) 



(37) 



The SBT being an involution, / = / so that //(r/„) = /(r/„). 
Much like the previous set of equations had introduced a discrete 
ki n grid, a discrete r&, grid is obtained for the radial component. 
Since Eqs. (34) and (36b) can be used to compute / and // for 
any value of r and k, they can in particular be used to compute 
f{ri n ) and fi(ri> n ) where V does not have to match /. The SBT and 
its inverse can then be expressed only in terms of series: 



In practical applications, for a given value of / only a limited 
number N max of fi{ki n ) and /(r/„) coefficients can be stored so 
that r lNmax = R and k lNmax = K h Since r in is defined by r ln = f^, 
for n - N max R and Ki are bound by the following relationship: 

qw = KiR. (39) 

Therefore, from the value of R one can easily determine the ap- 
propriate value for K\. Furthermore, because of the boundary 
conditions on // and /, the series (38a) and (38b) can now be 
truncated to N max terms. 

This truncation allows the use of a very convenient matrix 
formalism to represent the transform. In their explicit form, the 
two truncated sums can be rewritten as: 



Mk'n) = 
firvn) = 



Yj Kn P )Ki p ji{r lp k Vn \ (38a) f^) = Kf J f(r lp ) 



^/2n . / qi p qi>n 



2 iww/W- (38b) f{rvn ) = R-i y fi(ki P )-^^-ji l q -^) . 

p=l pTl Jl+Mlp) 



(40) 



(41) 



With these equations one can easily compute the SBT and its From these equations, a transform matrix T 11 ' can be defined as: 
inverse without the need of evaluating any integral. Furthermore 



only discrete values of / and / respectively sampled on r\ n and 
ki n are required. 



V2tt .Ai> P qiq. 
-Ji(- ) 



77+1 (qi q ) qiNnu 



(42) 



pq 
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This matrix allows the computation of // on any grid ky n from 
the values of / sampled on r\ n \ 



fi(kn) 
fi(kn) 

M k l'N max ) 



f(rn) 
f(ra) 

Rri Nm J 



(43) 



Reciprocally, the inverse the values of / can be computed on any 
?>„ grid from fi sampled on ky n using the exact same matrix: 



f{rn) 
f(rn) 



fl 3 



Mil) 
Mkw) 



(44) 



The discrete SBT is defined by the set of equations (43) and 
(44). 

A simplified form of the transform could have been defined 
only for / = V so that T lv = T l . However, keeping the distinction 
between the order of the transform / and the order of the grid on 
which the results are provided V will be crucial to the implemen- 
tation of the DSFBT. Indeed, the order of the grid on which the 
function is sampled has to match the order of the transform but 
the resulting transform coefficients do not. Therefore, it will be 
possible to compute the result of the inverse SBT of any order / 
on a grid of order Iq so that only one radial grid of order Iq will 
be required. Nevertheless, for the direct transform, if the field 
is sampled on the radial grid of order / , only the transform of 
order /q can be computed. An additional result is required to be 
able to relate the SBT of different orders between them. This is 
achieved by combining relations (33) and (34): 



flikm) = 



fir 

n Jo 
12 



fl (klom)Pl mjlo(klom r ) 



ji(ki n r)r dr, 



|2f . r R 2 

\ - ) ,fi (ki om )pi om ji (ki om r)ji(ki n r)r dr, 

v - 2 f 1 • 

ti JJ 0+ Mom) JO 

2 



= Yjflo(kl m)- 



jl (qi Q mX)jl(qinX)x dx, 



--1 "io 
ill' 



A+iteow) 



w l ° l 



(45) 



(46) 



where the weights W nm are defined as: 

W nl = I jlMkmX)jl{qinX)X 2 dx. 
JO 

The final expression in Equation 45 is an important result 
which shows that the SBT of a given order can be expressed as 
the sum of the coefficients obtained for a different order of the 
transform, with the appropriate weighting. This means we can 
convert the Spherical Bessel coefficients of order / into coeffi- 
cients of any other order /, which considerably speeds up calcu- 
lations for the SBT. It is also worth noticing that the weights W l n l m 
are simply geometric terms, i.e. independent of the field and can 
thus be tabulated. 

We note that this approach is an extension of the Discrete 
SBT presented in Lemoine (1994) using spherical Bessel func- 
tions and where the transform in Lemoine (1994) can be consid- 
ered as a special case where / = V . 



4.2. A Spherical 3D grid compatible with a Discrete 
Spherical Fourier-Bessel spectrum 

As presented in section 2.1, the SFB transform is the composi- 
tion of a SHT for the angular component and a SBT for the radial 
component. Since these two transforms can commute, they can 
be treated independently and by combining discrete algorithms 
for both transforms, one can build a DSFBT. In this work, the an- 
gular part of the transform is implemented using the HEALPix 
pixelisation scheme while the radial component uses the discrete 
SBT algorithm presented in Section 4. 1 . The choice of these two 
algorithms introduces a discretisation of the SFB coefficients as 
well as a 3D gridding of the density in physical space. 

The SFB coefficients fi m (k) are defined by Equation (1) for 
continuous values of k. Assuming a boundary condition on the 
density field /, the discrete SBT can be used to discretise the 
values of k. The Discrete SFB coefficients are therefore defined 
as: 



flmn — flmikin), 



(47) 



for < / < L max , -I < m < I and 1 < n < N max . These discrete 
coefficients are simply obtained by sampling the continuous co- 
efficients on the ki n grid introduced in the previous section. 

To this discretised SFB spectrum corresponds a dual grid of 
the 3D space defined by combining the HEALPix pixelisation 
scheme and the discrete SBT. 

Used for the computation of the SHT, the HEALPix scheme 
maps the sphere using curvilinear quadrilaterals of different 
shapes but equal area. Although other pixelisation schemes for 
the sphere such as IGLOO or GLESP could have been used, the 
choice of HEALPix is essentially motivated by the homogeneity 
of the pixelisation and by its comprehensive software package. 
For a given value of r, the field f(r, 6, (p) can therefore be sam- 
pled on a finite number of points using HEALPix. 

The radial component of the transform is conveniently per- 
formed using the discrete SBT. Indeed, this algorithm intro- 
duces a radial grid compatible with the discretised k\ n spectrum. 
Although this radial grid r/„ depends on the order / of the SBT, 
it will be justified in the next section that only one grid r/ oW is re- 
quired to sample the field along the radial dimension. The value 
of Iq is set to because in this case the properties of the zeros 
of the Bessel function ensure that r§ n will be regularly spaced 
between and R: 

r 0n = -^R. (48) 

For given values of 6 and 0, the field f(r, 0, (p) can now be sam- 
pled on discrete values of r = r§ n . 

Combining angular and radial grids, the 3D spherical grid is 
defined as a set of N max HEALPix maps equally spaced between 
and R. An illustration of this grid is provided on Fig. 3 where 
only one quarter of the space is represented for clarity. 

With this 3D grid, it will be possible to compute back and 
forth the SFB transform between a density field and its SFB co- 
efficients. The details of the actual algorithm are provided in the 
next section. 



4.3. Algorithm for the Discrete Spherical Fourier-Bessel 
Transform (DSFBT) 

Although the SHT and the SBT formally commute, some prac- 
tical considerations on the order between the two operations are 
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-0.2- 



-0.4- 




if 




Fig. 3. Representation of the spherical 3D grid for the DSFBT 
(R = 1 and N max = 4) as described in Section 4 where the axes 
represent x, y, z and the grid points correspond to the centre of 
Healpix pixels and the radial grid is defined in Section 4.1. 



to be taken into account for their actual implementation. Here, a 
detailed algorithm of both the direct and inverse DSFBT is pro- 
vided. 



4.3.1. Direct Transform 

Given a density field / sampled on the spherical 3D grid, the 
SFB coefficients fi mn are computed in three steps: 

1) For each n between 1 and N max the SHT of the HEALPix 
map of radius r/ oW is computed. This yields // m (r/ „) coeffi- 
cients. 

2) The next step is to compute the SBT of order /o from the 
fimifi^n) coefficients for every (/, m). This operation is a sim- 
ple matrix product: 









' fim(r k i) ' 


/V < / < L max 




flm(n 2) 


\y -I < m < I ' 




= ~w 


Jlm(n N max ). 



This operation yields f l lm (ki Qn ) coefficients which are not yet 
SFB coefficients because the order of the SBT Iq does not 
match the order of the Spherical Harmonics coefficients /. 
An additional step is required. 

The last step required to gain access to the SFB coefficients 
fi mn is to convert the Spherical Bessel coefficients for order 
/o to the correct order / that matches the Spherical Harmonics 
order. This is done by using relation (45): 



(49) 



3) 



(V < / < L max 
V -/ < m < I 
|V 1 < n < N max 



2W 



■hi 



flmikn) = Yftihp) .(50) 



p=\ 



where W n % are defined by Eq. (46). This operation finally 
yields the f imn = fi m (ki n ) coefficients. 

4.3.2. Inverse Transform 

Let fi mn be the coefficients of the SFB transform of the density 
field /, where the f mn coefficients can be calculated for a con- 
tinuous density field using the direct transform described above 
or for a galaxy or halo catalogue using the public code 3DEX 
(Leistedt et al. 2011). Note that the 3DEX code can account for 
masked regions of missing data. The reconstruction of / on the 
spherical 3D grid requires two steps: 

1) First, from the // m „, the inverse discrete SBT is computed 
for all / and m. Again, this transform can be easily evaluated 
using a matrix product: 





" fim(n i) ' 




flml 


< / < L m ax 


flm(n 2) 


Jill 


flml 


-I < m < I ' 




= W 






Jlm(n N max ). 




SlmNmax- 



(51) 

Here, it is worth noticing that the matrix T ll ° allows the eval- 
uation of the SBT of order / and provides the results on the 
grid of order Iq. 

2) From the spherical harmonics coefficients // m (r/ „) given 
at specific radial distances r/ „ it is possible to compute 
the inverse SHT. For each n between 1 and N ma x the 
HEALPix inverse SHT is performed on the set of coefficients 
{fim(ri Q n)}i,m- This yields N ma x HEALPix maps which consti- 
tute the sampling of the reconstructed density field on the 3D 
spherical grid. 

5. Wavelet decomposition of a test density field 

To illustrate the wavelet transform described in Section 3, a set 
of SFB Coefficients was extracted from a 3D density field using 
the DSFBT described in in Section 4. The test density field was 
provided by a cosmological N-body simulation which was car- 
ried out by the Virgo Supercomputing Consortium using com- 
puters based at Computing Centre of the Max-Planck Society in 
Garching and at the Edinburgh Parallel Computing Centre. We 
use their large box simulations 1 . 

We first compute the SFB coefficients of the test density field 
by sampling the Virgo density field on the 3D grid illustrated in 
Figure 3, for nside = 2048, / max = 1023 and n m?LX = 512. In 
order to perform the SFB decomposition, we place ourselves in 
the middle of the box, and calculate the SFB coefficients out to 
r = 479/2 h~ l Mpc, setting the over-density field to zero outside 
of this spherical volume. In order to better visualise the data, 
we present in Figure 4(a) the data in a cube spanning half the 



1 http://www.mpa-garching.mpg.de/Virgo/VLS.html, a 
ACDM simulation at z = 0, which was calculated using 512 3 par- 
ticles for the following cosmology : Q m = 0.3, Q A = 0.7, H = 



lOkms Y Mpc 
length 



, cr 8 = 0.9. The data cube provided is 479h Mpc in 
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size of the original box (i.e. 239.5/T 1 Mpc), where the field has 
been reconstructed from the SFB coefficients, using the discrete 
inverse transform we discuss in Section 4. 

From the SFB coefficients we calculate the wavelet de- 
composition, which yields the SFB coefficients of the various 
wavelet scales and the smoothed density field. Using the inverse 
DSFBT, the actual wavelet coefficients can be retrieved in the 
form of 3D density fields. These density fields corresponding to 
different wavelet scales are shown on Figure 4(b) to 4(e). The 
smoothed density field which arises from the wavelet decompo- 
sition is also shown in Figure 4(f), which is simply given by: 

CO = {a) - (b) - (c) - (d) - (e). (52) 




(a) Reconstructed density from (b) First wavelet scale 

the initial SFB coefficients 




(e) Fourth wavelet scale (f) Smoothed density 



Fig. 4. Spherical 3D wavelet decomposition of a density field for 
a box with side 239.5/T 1 Mpc. Panel (a) shows the reconstructed 
density field directly using the initial SFB coefficients. Panels 
(b)-(e) show the first four wavelet scales, while panel (f) shows 
the smoothed density field. 



6. Application to wavelet hard thresholding 
(denoising) 

In this section, a noise removal application based on wavelet 
thresholding is presented as an example of a potential use for 
the Isotropic Undecimated Spherical 3D Wavelet transform. 

6.1. Thresholding 

Many wavelet filtering methods have been proposed in the last 
twenty years. Hard thresholding consists of setting to all 
wavelet coefficients which have an absolute value lower than a 
threshold Tj (non-significant wavelet coefficient): 

# = ( W J* if \Wj, k \>Tj, 
h | otherwise, 

where is a wavelet coefficient at scale j and at spatial posi- 
tion k. 

Soft thresholding consists of replacing each wavelet coeffi- 
cient by the value w where 

f sgn(w jtk )(\ w jtk | -Tj) if | w# \> Tj, 
j,k 1 otherwise. 

This operation is generally written as: 

Wj, k = soft(w^) = sgn(w jtk )(\ w jJc I -Tj)+, (53) 

where (x)+ = max(0, x). 

In the case of Gaussian noise with a given standard deviation 
<r, a reasonable choice for the threshold Tj is Tj = Ko~j, where 
j is the scale of the wavelet coefficient, <x 7 is the noise standard 
deviation at the scale j, and K is a constant generally chosen 
between to 3 and 5. The standard deviation a j depends only on 
the chosen wavelet function and the noise level cr. More details 
can be found in Starck et al. (2010). 

Other threshold methods have been proposed, like the uni- 
versal threshold (Donoho & Johnstone 1994; Donoho 1993), or 
the SURE (Stein Unbiased Risk Estimate) method (Coifman & 
Donoho 1995), but they generally do not yield as good results 
as the hard thresholding method based on the significant coef- 
ficients. Concerning the threshold level, the universal threshold 
corresponds to a minimum risk. The larger the number of pixels, 
the larger the risk is, and it is normal that the threshold T depends 
on the number of pixels (T = o~j ^2 log ft, n being the number of 
pixels). The Kcr threshold corresponds to a false detection prob- 
ability, the probability to detect a coefficient as significant when 
it is due to the noise. The 3<x value corresponds to 0.27 % false 
detection. 

6.2. Denoising experiment 

We test the hard thresholding algorithm by using it to remove 
artificially added Gaussian noise on the Virgo density field (the 
simulation is described in Section 5). To conduct this experi- 
ment, the original density field taken from the n-body simula- 
tion has been used to compute an initial set of SFB coefficients 
(with l max = 1023 and N max = 512) and out to r = 479/T 1 Mpc. 
Figure 5(a) shows a slice of the 3D density field reconstructed 
from these coefficients. A Gaussian noise was then added to the 
SFB coefficients, the reconstruction of this noisy density field is 
shown on Figure 5(b). The hard thresholding algorithm was sub- 
sequently applied to the noisy SFB coefficients using 5 wavelet 
scales. The resulting density field is reconstructed on Figure 
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(c) Density field after wavelet hard thresholding 
Fig. 5. Wavelet Hard thresholding applied to a test density field 



5(c). The residuals are shown on Figure 5(d). The wavelet anal- 
ysis means we can successfully remove the noise we artificially 
added on entry , without much loss to the large scale structure, 
though some of the smaller structures are removed. 



7. Conclusion 

Modern cosmology requires the analysis of 3D fields on large 
areas of the sky, i.e. where the field is best viewed in spheri- 
cal coordinates. In this configuration, a Spherical Fourier Bessel 
(SFB) transform is the most natural way to statistically anal- 
yse the field. Wavelet transforms have been shown to be ideally 
suited for cosmological fields, which tend to be sparse in wavelet 



space. The wavelet transform can be used e.g. for denoising, but 
there is yet no 3D wavelet transform in spherical coordinates. 

We present in this paper a new 3D spherical wavelet trans- 
form, based on the undecimated wavelet transform (UWT) de- 
scribed in (Starck et al. 2006). In order to perform operations 
on the wavelet transforms (such as denoising), we require a dis- 
crete version of the SFB transform for both the direct and in- 
verse transforms. We show a novel way to perform such a fast 
Discrete Spherical Fourier-Bessel Transform (DSFBT) based on 
both a discrete Bessel transform and the HEALPIX angular pix- 
elisation scheme. 

Using the 3D wavelet transform and the DSFBT, both intro- 
duced in this paper, we denoise a test large scale structure data 
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set, taken from the Virgo large box simulations 2 . We find we can 
satisfactorily remove artificially added Gaussian noise without 
much loss to the large scale structure. All the algorithms pre- 
sented in this paper are available for download as a public code 
called MRS3D at http : // j starck . free . f r/mrs3d . html. 
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Spherical Fourier-Bessel decomposition and discretisation scheme. The 3D 
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Appendix A: Spherical Fourier-Bessel Transform 
and 3D Convolution Products 

A.1. Spherical Fourier-Bessel Transform and relation to the 
3D Fourier Transform 

In order to have a better understanding of the SFB coefficients 
and of how to use them to perform filtering, the SFB transform 
can be related to the 3D Fourier transform. We follow a similar 
definition as the one presented in Baddour (2010), but using our 
conventions for the different transforms. The following conven- 
tion will be used for the Fourier transform: 



F(k) 



2tt) 3 / 



f(r)e~ ikr dr, 



(A.1) 



where F denotes the Fourier transform of /. This formulation 
does not assume any coordinate system. However, to relate this 
transform to the SFB transform, it is possible to express this 
equation in spherical coordinates using the following expansion 
for the Fourier kernel: 



■i) l Mkr)Y?(6 r9 <f>r)Y™(6 k , (A.2) 



1=0 / 



where (k, 6k, <f>k) an d (r, 6 r , r ) are respectively the spherical co- 
ordinates of vectors k and r. 

Substituting this expression for the kernel in the definition of 
the 3D Fourier transform yields: 



F(k,e k ,<p k ) = —= VV(-/)'/(^r,^ 



) 



X Mkr)Y™(9 r , <p r )Y?(6 k , (f> k )dQr l dr, 
X MkrjYjW, 4>r)d£lr 2 dr] Yf(6 k , </> k ), 

GO I 

[(-i) l fUk)]Y?(e k ,<p k ). 

1=0 m =-l 



(A3) 



http : //www . mpa- garching . mpg . de/Virgo/VLS . html 



In the last equation, the expression of the Spherical 
Harmonics Expansion of F(k, 6%, (pk) for a given value of k can be 
recognised from Eq. (4b). In the Fourier space, the (-i) l fi m (k) are 
the Spherical Harmonics coefficients of F on a sphere of given 
radius k. In other words, the Spherical Harmonics coefficients 
Fi m (k) of the 3D Fourier transform F(k, 6k, (pk) on a sphere of 
given radius k in Fourier space are the SFB coefficients fi m (k) 
for the same value k but multiplied by factor {-if. 

The relationship between 3D Fourier transform and SFB 
transform is therefore very simple. The SFB transform can be 
sought of as a mere Fourier transform in spherical coordinates. 
In the next sections, this relationship will be used to derive con- 
volution and filtering relations for the SFB transform using the 
well known relations verified by the Fourier transform. 
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A.2. Expression of a 3D convolution product using the SFB 
Transform 

A prerequisite to the establishment of filtering relations is the 
expression of a 3D convolution product in terms of SFB coef- 
ficients. Let v(r, r , <p r ) be the 3D convolution of f(r, 6 r , <p r ) and 
u(r, 6 r , (p r ). Then the 3D Fourier transform of v verifies: 

V(k,6 k ,(f> k ) = T{f*uKk,6 k ,<f> k ), 



= ^(2^F(k,e k ,(/> k )U(k,e k ,(f> k ). (A.4) 

where f denotes the 3D Fourier transform. From Eq. (A. 3) the 
expression of the 3D Fourier transform in spherical coordinates 
is known in terms of SFB coefficients. Applying this relationship 
to V(k, 6 k , (p k ) in the last equation yields: 



(-i) l vim(k) = f f ^2^F(k,6 k ,(f> k )U(Ke k ,(f> k ) 
Jo Jo 



x Y™(6 k ,(f> k )sm((f> k )d(f> k d6 k . (A.5) 
Then, by applying Eq. (A.3) to F and U one gets: 



v lm (k) = (i) l ^^(2n? flZZ (-rffi'mWYfiOk, 

J J /'=0m'=-/' 

oo I" 

x J] J] (-O r W(*)<(0 t ,&) 



x Y™(e k ,</> k )sm(<t> k )d<f> k de k , 

oo /' 

= (o'V(2^Z Z (-o7pm-(*) 



l'=0 m'=-l' 

I" 

X 

Z"=0 m"=-V 



X 



Z Z (- i y" U l"m"(k) 

Z"=0 m"=-l" 

Jj if (ft, (f> k )Yf(6 k , <f> k )Y»(fi k , (f> k )dQ k . (A.6) 



The last integral over the two angular variables can be expressed 
as a Slater integral (which is a special case of the Gaunt integral) 
defined as: 

c l '\l,mj',m') = Jj YfiQ, $)Yf' (0, (p)Yfr m> (0, (A.7) 

The Slater integrals are only nonzero for \l - l'\ < I" < I + V 
which simplifies the expression of vi m (k). 

The SFB transform of the 3D convolution product is there- 
fore: 



(f*u) lm (k) = (o'Vc^Z Z (-o7^(*) 

Z'=0 m'=-Z' 

i+r 

Z c /,, (/,m,/ , ,m , )(-0 /,, ^m-m'W. (A.8) 



i+r 

x 

i"=\i-v\ 
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